Overview

Plant part and environmental context shape symbiotic microbial compositions, yet the interaction of these two variables are seldom considered. We analyzed epiphytic microbes from nine Hibiscus tiliaceus (Hawaiian: hau) trees across a steep environmental gradient within a single Hawaiian watershed in Waimea Valley on the north shore of O‘ahu, Hawai‘i. At each location we sampled 8 microhabitats: leaves, petioles, axils, stems, roots, and litter from the plant, as well as surrounding air and soil.

The composition of microbial communities is driven primarily by microhabitat (i.e., above vs. below ground communities), this variable predicted more than twice the compositional variance in bacteria compared to fungi. Fungal community compositions, in contrast, were more responsive to gradient dynamics than bacteria, and there were differences between spatial dynamics of fungal communities associated with different microhabitats that correlate with distribution patterns and range size. Within plants, microbes were compositionally nested with aboveground communities containing a subset of diversity found belowground. Our findings suggest potential differences in the underlying mechanisms shaping communities of fungi and bacteria associated with plants, and indicate an interaction between assembly mechanisms working simultaneously on different spatial scales.

if (!require('knitr')) install.packages('knitr'); library('knitr')
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.align='center')

# Load in packages
if (!require("pacman")) install.packages("pacman"); library(pacman) # for rapid install if not in library
if(!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager") #Load packages from Bioconductor

pacman::p_load('FedData', 'latticeExtra', 'scales', 'viridis', 'ggmap', 'phyloseq', 'raster', 'rgdal', 'RColorBrewer', 'ggplot2', 'vegan', 'ggplot2', 'phyloseq', 'bipartite', 'plotrix', 'viridis', 'lattice', 'fossil', 'plyr', 'devtools', 'ggpubr', 'gridExtra', 'cowplot')

devtools::install_bitbucket("graumannlabtools/multipanelfigure")

Site Map

The site where the study was conducted: Waimea Valley on the north shore of O‘ahu, Hawai‘i.

# Code by Austin Greene
# Edited by Chris Wall

# Purpose: Generate a map of Waimea Valley, Oahu with elevation data, 
# and site locations which are scaled the mean-annual precipitation 

# Load in physeq object
physeq1=readRDS("data/physeq1.gz")

# Getting base map via Stamen of location within bounding box. 
# Using Stamen maps which do not require a Google API key
map3 <- get_map(location = c(left = -158.070631, bottom = 21.598356, right = -157.998464, top = 21.655101), zoom=13, source="stamen", maptype = c("terrain-background"))
Waimea_map_3 <- ggmap(map3) # Turn into ggmap graphical object

# Generating extents of map to pull elevation
bb <- attr(map3, "bb") # Mapp attribute object, from which we extract map extents 
extentB <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon, bb$ll.lat, bb$ur.lat),
                               proj4string = "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")

# Confirm bounding box for our elevation data
# ggmap(map3) + geom_polygon(data=extentB, fill=NA, aes(x=long,y=lat), color="red", size=3)

# Download elevation data tiles
elev_waimea <- get_ned(template = extentB, label = "ned_waimea", res="1", force.redo = F)

# Make df of elevation data with lat and long included
elev_raster_df <- raster::as.data.frame(elev_waimea, xy=TRUE) #xy True 

# Make a dataframe of geographic location and mean annual precip at each site (from existing metadata)
geo_p=cbind(sample_data(physeq1)$Long, sample_data(physeq1)$Lat, sample_data(physeq1)$rain, sample_data(physeq1)$FieldSite)
geo3 = data.frame(geo_p) #we convert it to a dataframe
colnames(geo3)=c("Long", "Lat","Rain", "Site") # Repair column names
geo3 <- subset(geo3, Site != 9) # Remove Site 9 
geo3$Site <- factor(geo3$Site) # Convert sites to factors, not critical

# Plot elevational data and site locations scaled by precipitation on top of existing ggmap of Waimea Valley
Waimea_map_6_sizescaled <- Waimea_map_3 + # Existing map
  geom_raster(data = elev_raster_df, aes(x=x, y=y, fill=elev_raster_df$ned_waimea_NED_1)) + # Raster elevation data
  geom_point(data = geo3, aes(x=Long, y=Lat, size=Rain), fill="white", color="black", pch=21, alpha=0.5) + # Points for sites, scaled by mean-annual precipitation
  coord_cartesian() +  # Set to cartesian coordinates
  scale_fill_viridis(option = "plasma", alpha = 0.7) + # Set fill colors to be colorblind-friendly
  coord_fixed(1.3) + # Aspect ratio of cartesian coordinates is 1.3
  labs(x="Longitude", y="Latitude", size="Precipitation (mm/year)", fill="Elevation (m)") + # Labels
  theme_classic() # Classic minimalist theme

Waimea_map_6_sizescaled # Print the plot

Figure 1. Site map with sampling locations

Distance decay: Bacteria and fungi

Inspect the bacterial and fungal communities using rarified and square-root transformed counts un Bray-Curtis dissimilarity matrix. There is a significant correlation between the 2 communities, with each showing similar patterns of dissimilarity and similarity within individual samples.

# Code by Anthony
# modifed by Chris Wall
# To plot mantel correlation between fungi and Bacteria

#read normalized Bacteria Phylloseq object
PS_normal=readRDS("data/PS_normal.gz")
PS_normal<- subset_taxa(PS_normal, c(Kingdom != "unknown", Kingdom != "Eukaryota"))

#read normalized Fungi Phylloseq object
fungal_PS_normal=readRDS("data/fungal_PS_normal.gz")

#sqrttransform
sqrtfun=transform_sample_counts(fungal_PS_normal, function(x) x^0.5)

#remove the sample missing from Bacteria
sqrtfun=subset_samples(fungal_PS_normal, CollectionID !="WMEA_01223_Pl_1")

#calculate bray-curtis distance 
funhel.dist=phyloseq::distance(sqrtfun, "bray") 

#remove this sample wich was low abundance in fungi
PS_normal=subset_samples(PS_normal, CollectionID !="WMEA_01223_Pl_1")

#sqrttransform
sqrt=transform_sample_counts(PS_normal, function(x) x^0.5)

#calculate bray curtis distance
hel.dist=phyloseq::distance(sqrt, "bray") 

#Calculate mantel
m=mantel(hel.dist,funhel.dist)

#make a dataframe for ggplot
mantel=as.data.frame(cbind(c(hel.dist), c(funhel.dist)))

# make plot
ggplot(mantel,aes(mantel$V1,mantel$V2))+
  geom_point()+
  geom_smooth()+
  theme_bw()+
  xlab("Fungal Dissimilarity")+
  ylab("Bacteria Dissimilarity")+
  annotate("text", x=.4, y=.3, label= "r=0.434"~italic("p")~"=0.001", size=3)+
  theme(plot.title = element_text(hjust = 1))
Figure. Bacterial and Fungal community mantel tests.

Figure. Bacterial and Fungal community mantel tests.

dev.print(pdf, "figures/FigS3. Community.mantel.pdf", width=4, height=4)
dev.off()

Abundance-occupancy

Generate abundance-occupancy graphs with bacteria and fungal data with site locations (habitats), which are scaled the mean-annual precipitation. Now generate plots for Bacteria and Fungi:
1. Mean ASV range for ASVs or Genus levels
2. Number of Habitats ASV is Present in vs ASV range
3. Mean abundance per sample (when present) occurrence vs. ASV range

Statistical tests in each figure show the significance of relationships using correlations tests and fitting line from lsfit.

# Code by Anthony; modified by Jared Bernard and Feresa Corazon
# Purpose: Generate abundance occupancy graphs with bacteria and fungal data and site locations which are scaled the mean-annual precipitation 

# Function
phyloseq_standardize_otu_abundance <- function(physeq, method="total", ...){  
  trows <- phyloseq::taxa_are_rows(physeq)  ## Check the orientation of the OTU table
  if(trows == TRUE){ marg <- 2 } else { marg <- 1 }
  comm <- as(object = phyloseq::otu_table(physeq),
             Class = "matrix")  ## Extact OTU table
  comm_std <- vegan::decostand(comm, method,
                               MARGIN = marg, ...)  ## Standardize community table
  phyloseq::otu_table(physeq) <- phyloseq::otu_table(comm_std,
                                                     taxa_are_rows = trows) ## Replace old otu_table with the new one
  return(physeq)
}


#Section 1A. Calling all data and values for bacteria
##################################
## Bacterial Abundance Occupancy ##
##################################

rarPhySeq=readRDS("data/rarPhySeq.gz")
rarPhySeq <- subset_taxa(rarPhySeq, Kingdom != "unknown")   #Remove unclassified kingdoms
rarPhySeq <- subset_taxa(rarPhySeq, Kingdom != "Eukaryota") #Remove eukaryotes

#Calculate distance to shore
sample_data(rarPhySeq)$Shore_dist = pointDistance(cbind
    (sample_data(rarPhySeq)$Long,sample_data(rarPhySeq)$Lat),
    c(-158.062848, 21.640741),lonlat = TRUE)

# Merge by sample type
agg=merge_samples(rarPhySeq, "SampleType")

#How many habitats is each ASV found?
#Standardize to convert to presence absence
habitats = colSums(otu_table(phyloseq_standardize_otu_abundance(agg,
    method = "pa")))

#Merge by site location
sites = merge_samples(rarPhySeq, "FieldSite")

#Standardize to convert to presence absence and calculate site total
site = colSums(otu_table(phyloseq_standardize_otu_abundance(sites,
    method = "pa")))

#Convert the otu table to presence absence in order to calculate range size
binarysite = phyloseq_standardize_otu_abundance(sites, method = "pa")

#Multiply each ASV by its distance to shore
range = otu_table(binarysite)*sample_data(binarysite)$Shore_dist

#Convert 0 to NA
is.na(range) <- range == 0

#Calculate min and max distance to shore
rangespan = apply(range,2,range, na.rm=TRUE)

#Subtract min distance from max distance
rangespan = rangespan[2,]-rangespan[1,]

#Do same thing but log transform range
nozerorangespan = rangespan
is.na(nozerorangespan) <- nozerorangespan == 0

#Calculate total abundance
abund = colSums(otu_table(rarPhySeq))     

#Calculate how many samples an ASV is present
occupancy = colSums(otu_table(phyloseq_standardize_otu_abundance(rarPhySeq,
    method = "pa")))

#Calculate mean abundance per sample (where present)
meanabund = otu_table(rarPhySeq)
is.na(meanabund) <- meanabund == 0
meanabund = colMeans(meanabund, na.rm = TRUE)

#Calculate per site habitat diversity (where present)
habpersite = cbind(colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s1"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s2"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s3"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s4"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s5"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s6"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s7"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s8"), method = "pa"))),
                 colSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeq, FieldSite=="s10"), method = "pa"))))
is.na(habpersite) <- habpersite == 0
habpersite = rowMeans(habpersite, na.rm = TRUE)

##################################
## Bacterial Abundance Occupancy at Genus Level ##
##################################

genus_sites = tax_glom(sites, taxrank = "Genus")

#Standardize to convert to presence absence and calculate site total
genus_site = colSums(otu_table(phyloseq_standardize_otu_abundance(genus_sites,
    method = "pa")))

#Convert the ASV table to presence absence in order to calculate range size
genus_binarysite = phyloseq_standardize_otu_abundance(genus_sites,
    method = "pa")

#Multiply each ASV by its distance to shore
genus_range = otu_table(genus_binarysite)*
    sample_data(genus_binarysite)$Shore_dist

#Convert 0 to NA
is.na(genus_range) <- genus_range == 0

#Calculate min and max distance to shore
genus_rangespan = apply(genus_range,2,range, na.rm = TRUE)

#Subtract min distance from max distance
genus_rangespan = genus_rangespan[2,]-genus_rangespan[1,]

#Section 1B. Calling all data and values for fungi
###############################
## Fungal Abundance Occupancy##
###############################

rarPhySeqF=readRDS("data/fungal_rarPhySeq.gz")

sample_data(rarPhySeqF)$Shore_dist = pointDistance(cbind
    (sample_data(rarPhySeqF)$Long,sample_data(rarPhySeqF)$Lat),
    c(-158.063625, 21.640741), lonlat = TRUE)

aggF = merge_samples(rarPhySeqF, "SampleType")

#How many habitats is each ASV found?
#Standardize to convert to presence absence
habitatsF = colSums(otu_table(phyloseq_standardize_otu_abundance(aggF,
    method = "pa")))

#Merge by site location
sitesF = merge_samples(rarPhySeqF, "FieldSite")

#Standardize to convert to presence absence and calculate site total
siteF = colSums(otu_table(phyloseq_standardize_otu_abundance(sitesF,
    method = "pa")))

#Convert the ASV table to presence absence in order to calculate range size
binarysiteF = phyloseq_standardize_otu_abundance(sitesF, method = "pa")

#Multiply each ASV by its distance to shore
rangeF=otu_table(binarysiteF)*sample_data(binarysiteF)$Shore_dist

#Convert 0 to NA
is.na(rangeF) <- rangeF == 0

#Calculate min and max distance to shore
rangespanF = apply(rangeF,2,range.default, na.rm = TRUE)

#Subtract min distance from max distance
rangespanF = rangespanF[2,]-rangespanF[1,]

#Do same thing but log transform range
nozerorangespanF = rangespanF
is.na(nozerorangespanF) <- nozerorangespanF == 0

#Calculate total abundance
abundF = colSums(otu_table(rarPhySeqF)) 

#Calculate how many samples an ASV is present
occupancyF = colSums(otu_table(phyloseq_standardize_otu_abundance(rarPhySeqF,
    method = "pa")))

#Calculate mean abundance per sample (where present)
meanabundF = otu_table(rarPhySeqF)
is.na(meanabundF) <- meanabundF == 0
meanabundF = rowMeans(meanabundF, na.rm = TRUE)

#Calculate per site habitat diversity (where present)
habpersiteF = cbind(rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s1"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s2"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s3"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s4"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s5"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s6"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s7"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s8"), method = "pa"))),
                  rowSums(otu_table(phyloseq_standardize_otu_abundance(subset_samples(rarPhySeqF, FieldSite=="s10"), method = "pa"))))
is.na(habpersiteF) <- habpersiteF == 0
habpersiteF = rowMeans(habpersiteF, na.rm = TRUE)

##################################
## Fungal Abundance Occupancy at Genus Level##
##################################

genus_sitesF = tax_glom(sitesF, taxrank = "Genus")

#Standardize to convert to presence absence and calculate site total
genus_siteF = colSums(otu_table(phyloseq_standardize_otu_abundance(genus_sitesF,
    method = "pa")))

#Convert the ASV table to presence absence in order to calculate range size
genus_binarysiteF = phyloseq_standardize_otu_abundance(genus_sitesF,
    method = "pa")

#Multiply each ASV by its distance to shore
genus_rangeF = otu_table(genus_binarysiteF)*
    sample_data(genus_binarysiteF)$Shore_dist

#Convert 0 to NA
is.na(genus_rangeF) <- genus_rangeF == 0

#Calculate min and max distance to shore
genus_rangespanF = apply(genus_rangeF,2,range, na.rm = TRUE)

#Subtract min distance from max distance
genus_rangespanF = genus_rangespanF[2,] - genus_rangespanF[1,]

#Section 2. Making all plots 
######################################################
########## ALL PLOTS AND STATISTICAL TESTS ###########
######################################################

par(mfrow=c(3,2), mar=c(4,5,0.7,0.2))

#####################
#First plot: Bacterial rangespan (PLOT 1, row 1)
#Combine Genus and ASV tables
ranges <- data.frame(Genus=c(rep("Genus", times=length(genus_rangespan)), 
                        rep("ASV", times=length(rangespan))), 
                 combined_rangespan=c(genus_rangespan, rangespan))

#Set margin size for boxplot
# bactrange
boxplot(combined_rangespan~Genus, horizontal = TRUE, xaxt="n", yaxt="n", medlwd=1.5, lwd=0.8,
    vertical = TRUE, col = "gray48", frame = F, cex.lab = 0.6, ylab="Levels",
    main= "Bacteria", cex.main=0.8,
    xlab="", cex.axis = 0.6, data = ranges)
axis(side=1,labels=F, tck=-0.06)
axis(side=2,labels=c("ASV", "Genus"), cex.axis = 0.6, at=seq(2), tck=-0.06)
mtext(expression(bold('a')), side=3, line=-0.75, at=0, cex=0.5)

#####################
#Second Plot: Fungal rangespan (PLOT 2, row 1)
rangesF <- data.frame(Genus=c(rep("Genus", times=length(genus_rangespanF)), 
                             rep("ASV", times=length(rangespanF))), 
                     combined_rangespanF=c(genus_rangespanF, rangespanF))

#Set margin size for boxplot
# fungrange
boxplot(combined_rangespanF~Genus, horizontal = TRUE,  xaxt="n", yaxt="n", cex=0.4, medlwd=1.5, main="Fungi", cex.main=0.8, pch=16,
    col = "gray84", frame = F, cex.axis=0.8, ylab="", xlab="", lwd=0.8, data = rangesF)
axis(side=1,labels=F, tck=-0.06)
axis(side=2,labels=F, at=seq(2), tck=-0.06)
mtext(expression(bold('b')), side=3, line=-0.75, at=0, cex=0.5)

#compare ranges between bacteria and fungi genera
t.test(genus_rangespan,genus_rangespanF)

# Welch Two Sample t-test
# 
# data:  genus_rangespan and genus_rangespanF
# t = 5.7294, df = 869.79, p-value = 1.389e-08
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   467.2236 954.1336
# sample estimates:
#   mean of x mean of y 
# 3757.528  3046.850

#compare ranges between bacteria and fungi ASVs
t.test(rangespan, rangespanF)

# Welch Two Sample t-test
# 
# data:  rangespan and rangespanF
# t = 32.062, df = 13235, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   852.2800 963.2758
# sample estimates:
#   mean of x mean of y 
# 2077.959  1170.181

#####################
##Plot Bacteria: # of Habitats vs. ASV Range (PLOT 3, row 2)
# bact_habrange
plot(jitter(rangespan,20),jitter(habitats,3), xaxt="n", yaxt="n", ylim=c(0,10), pch=16,
    col="gray48", cex=0.2, xlab="", cex.lab=0.6,
    ylab="Number of Habitats \nASV is Present")+
    abline(lsfit(rangespan,habitats), col="black", lwd=1)
axis(side=1,labels=F, tck=-0.06)
axis(side=2, labels=T, tck=-0.06, cex.axis=0.6, at=seq(2,10, by=3))
mtext(expression(bold('c')), side=3, line=-0.75, at=0, cex=0.5)

#Calculate correlation
cor.test(rangespan, habitats)
#text(temp, expression(paste( italic("R^2 = 0.262 
#   p < 0.001"))))

#####################
#Plot Fungi: # of Habitats vs. ASV Range (PLOT 4, row 2)
# fung_habrange
plot(jitter(rangespanF,20),jitter(habitatsF,3), yaxt="n", xaxt="n", ylim=c(0,10), pch=16,
    col="gray74", cex=0.2, xlab="",  cex.axis=0.6, ylab="",
    abline(lsfit(rangespanF,habitatsF), col="black", lwd=1))
    axis(side=1,labels=F, tck=-0.06)
    axis(side=2,labels=F, tck=-0.06, at=seq(2,10, by=3))
mtext(expression(bold('d')), side=3, line=-0.75, at=0, cex=0.5)
    

#Calculate correlation
cor.test(rangespanF, habitatsF)
#text(temp, expression(paste( italic("R^2 = 0.451
#   p < 0.001"))))

#####################
# Plot Bacteria: Mean Abundance vs. ASV Range (PLOT 5,  row 3)
plot(jitter(rangespan,20),jitter(log(meanabund),3),col="gray48", cex=0.2, pch=16,
     xaxt="n", yaxt="n", ylim=c(0,10),
     cex.lab=0.6, xlab="ASV range (m)", ylab="Mean Abundance/Sample \n(Where Present)") +
axis(side=1,labels=T, tck=-0.06, cex.axis=0.6)
axis(side=2,labels=T, tck=-0.06, cex.axis=0.6, at=seq(2,10, by=3))
mtext(expression(bold('e')), side=3, line=-0.75, at=0, cex=0.5)

#Calculate correlation
cor.test(rangespan, meanabund)
#text(temp,expression(paste( italic("R^2 < 0.001
#   p = 0.554"))))

#####################
# Plot Fungi: Mean abundance vs. ASV Range (PLOT 6, row 3)
plot(jitter(rangespanF,20),jitter(log(meanabundF),3),col="gray74", cex=0.2, xaxt="n", ylim=c(0,10),
     yaxt="n", cex.lab=0.6, xlab="ASV range (m)", ylab="", pch=16) + 
  abline(lsfit(rangespanF,log(meanabundF)), col="black", lwd=1)
axis(side=1,labels=T, tck=-0.06, cex.axis=0.6)
axis(side=2,labels=F, tck=-0.06, at=seq(2,10, by=3))
mtext(expression(bold('f')), side=3, line=-0.75, at=0, cex=0.5)

#Calculate correlation
cor.test(rangespanF, meanabundF)
#text(temp,expression(paste( italic("R^2 = 0.012
#   p < 0.001"))))


dev.print(pdf, "figures/Fig4. ASVrangespan.pdf", width=6, height=5)
dev.off()

Additional taxonomy summary for abundance-occupancy data, used in Supplementary Figures 4 and 5

### TaxonomySummary ####
bacttax=tax_table(binarysite)
bactaxrange=data.frame(bacttax,rangespan,habitats)
#head(bactaxrange)

bacphybymedrange <- with(bactaxrange, reorder(Phylum, rangespan, median))
bacphybymedhab <- with(bactaxrange, reorder(Phylum, habitats, median))

### Bacteria phyla plot
par(mfrow=c(2,1), mar=c(4,5,1,1)+0.1)

# Bacteria rangespan
boxplot(bactaxrange$rangespan~bacphybymedrange, las=2, horizontal=TRUE, frame=F, cex.main=1,
        #main="Distribution of Bacteria Phyla", 
        col="gray48", cex=0.4, medlwd=1.5, lwd=1,
        xaxt="n", tck=-0.03, xlab="Range (m)", ylab="", cex.axis=0.3, cex.lab=0.8)
axis(side=1, labels=T, seq(0,5500, by=1000), cex.axis=0.6, cex.lab=0.6, tck=-0.02)
mtext(expression(bold('a')), side=3, line=-0.75, at=0, cex=0.8)

# Bacteria habitats
boxplot(bactaxrange$habitats~bacphybymedhab, las=2, horizontal=TRUE, frame=F, xaxt="n", tck=-0.03,
        col="gray48", xlab="Microhabitats", cex.axis=0.4, cex.lab=0.8, ylab="", cex=0.3, medlwd=1.5, lwd=1)
axis(side=1, labels=T, seq(0,8, by=1), cex.axis=0.6, cex.lab=0.6, tck=-0.02)
mtext(expression(bold('b')), side=3, line=-0.75, at=1, cex=0.8)

dev.copy(pdf, "figures/FigS4. BactTax.range.hab.pdf", width=6, height=8.5)
dev.off()


######### Fungi 
funtax=tax_table(binarysiteF)
funtaxrange=data.frame(funtax,rangespanF,habitatsF)
#head(funtaxrange)

funphybymedrange <- with(funtaxrange, reorder(Class, rangespanF, median))
funphybymedhab <- with(funtaxrange, reorder(Class, habitatsF, median))

### Fungi class plot
par(mfrow=c(2,1), mar=c(4,5,1,1)+0.1)

# Fungi rangespan
boxplot(funtaxrange$rangespanF~funphybymedrange,las=2, horizontal=TRUE, frame=F, cex.main=1,
        #main="Distribution of Fungal Classes", 
        col="gray84", cex=0.4, medlwd=1.5, lwd=1,
        xaxt="n", tck=-0.03, xlab="Range (m)", ylab="", cex.axis=0.4, cex.lab=0.8)
axis(side=1, labels=T, seq(0,5500, by=1000), cex.axis=0.6, cex.lab=0.6, tck=-0.02)
mtext(expression(bold('a')), side=3, line=-0.75, at=0, cex=0.8)

# Fungi habitats
boxplot(funtaxrange$habitatsF~funphybymedhab, las=2, horizontal=TRUE, frame=F, xaxt="n", tck=-0.03,
        col="gray84", xlab="Microhabitats", cex.axis=0.4, cex.lab=0.8, ylab="", cex=0.4, medlwd=1.5, lwd=1)
axis(side=1, labels=T, seq(0,8, by=1), cex.axis=0.6, cex.lab=0.6, tck=-0.02)
mtext(expression(bold('b')), side=3, line=-0.75, at=1, cex=0.8)

dev.copy(pdf, "figures/FigS5. FungTax.range.hab.pdf", width=6, height=8.5)
dev.off()

Nestedness plots

What is the nestedness for each community within the microhabitats of plant-parts and do they relate to each other among habitats (the sites of collection)?

#Code by Anthony
#To produce nestedness plots

#read in rarefied Bacteria data
rarPhySeq=readRDS("data/rarPhySeq.gz")
rarPhySeq<- subset_taxa(rarPhySeq, c(Kingdom != "unknown", Kingdom != "Eukaryota"))

#Merge samples by sample type
agg=merge_samples(rarPhySeq, "SampleType")

#convert to dataframe
aggtab=as.data.frame(otu_table(agg))

#calculate nested temperature
nested=nestedtemp(aggtab)

#read in Fungal data (not rarefied)
fungalphyseq1=readRDS("data/fungal_physeq1.gz")

#rarefy to 20000 sequences (drops 3 samples :())
funrarPhySeq=rarefy_even_depth(fungalphyseq1, sample.size=20000)

#merge by sample type
funagg=merge_samples(funrarPhySeq, "SampleType")

#convert to dataframe
funaggtab=as.data.frame(otu_table(funagg))

#calculate nested temp
funnested=nestedtemp(funaggtab)

#make a two row composite figure, reguce spacing between inner and outer margins
par(mfrow=c(2,1), mai=c(.1,1.,1,1), oma=c(.1,1,.1,1))
#plot the Bacteria nested figure, suppress taxon names
plot(nested, kind="incid", names=c(TRUE,FALSE), col=c("white", "gray48"), lwd=2, main="Bacteria", cex.main=0.8, cex.axis=0.8)
#add stats are margin text
mtext(expression(paste("Nested Temp=39.8, ", italic("p") ,"=0.001")), cex=0.7,)
mtext(expression(bold('a')), side=3, line=0, at=0, cex=0.8,)

#plot the fungal nested figure, suppress taxon names
plot(funnested, kind="incid", names=c(TRUE,FALSE), col=c("white", "gray84"), lwd=2, main="Fungi", cex.main=0.8, cex.axis=0.8)
#add stats as margin text
mtext(expression(paste("Nested Temp=42.5, ", italic("p") ,"=0.001")), cex=0.7)
mtext(expression(bold('b')), side=3, line=0, at=0, cex=0.8)

dev.print(pdf, "figures/Fig3. Nestedplots.pdf", width=6, height=6)
dev.off()

Mantel Tests

Tests here show the relationships between plant parts (or plant organs) and bacterial communities fungal communities across geographic distance. The exported table uses Hellinger square-root trasnformed abundance data in Mantel tests with p-values and Bonferroni corrected p-values.

## Mantel Tests and Slope Plots By Plant Part
## Chad Wilhite 04/25/19
## edited Chris Wall

####################################
#### Mantel Tests by SampleType ####
####################################

#
##
###
#### Bacteria Analysis

# Load RDS file of normalized bacterial reads
PS_normal = readRDS("data/PS_normal.gz") #"PS_normal"
PS_normal<- subset_taxa(PS_normal, c(Kingdom != "unknown", Kingdom != "Eukaryota"))

# Hellinger (square root) transform the data!
PS_normal.hell = transform_sample_counts(PS_normal, function(x) x^0.5)

# View SampleType names
# unique(get_variable(PS_normal.hell, sample_variables(PS_normal.hell)[18]))

### Subset out bacterial data by SampleType

# Subset bacterial data by plant part - the hard way... xD
stem_otu = subset_samples(PS_normal.hell, SampleType=="Stem")
root_otu = subset_samples(PS_normal.hell, SampleType=="Root")
air_otu  = subset_samples(PS_normal.hell, SampleType=="Air")
leaf_otu = subset_samples(PS_normal.hell, SampleType=="Leaf")
soil_otu = subset_samples(PS_normal.hell, SampleType=="Soil")
litt_otu = subset_samples(PS_normal.hell, SampleType=="Litter")
axil_otu = subset_samples(PS_normal.hell, SampleType=="Axil")
peti_otu = subset_samples(PS_normal.hell, SampleType=="Petiole")


### Generate Dissimilarity Matrices for bacterial OTU data

# Generate a dissimilarity matrix for grouped and individual plant part bacterial OTU data 
# adding "phyloseq::" to specify to run this through phyloseq and not other packages.

all_otu.dist  = phyloseq::distance(PS_normal.hell, "bray")
stem_otu.dist = phyloseq::distance(stem_otu, "bray") 
root_otu.dist = phyloseq::distance(root_otu, "bray") 
air_otu.dist  = phyloseq::distance(air_otu,  "bray") 
leaf_otu.dist = phyloseq::distance(leaf_otu, "bray") 
soil_otu.dist = phyloseq::distance(soil_otu, "bray")
litt_otu.dist = phyloseq::distance(litt_otu, "bray")
axil_otu.dist = phyloseq::distance(axil_otu, "bray") 
peti_otu.dist = phyloseq::distance(peti_otu, "bray")


### Create correctly ordered sized geographic distance dissimilarity matrix

#Gather bacterial spatial data from the phyloseq object
all.geo  = cbind(sample_data(PS_normal.hell)$Lon, sample_data(PS_normal.hell)$Lat)
stem.geo = cbind(sample_data(stem_otu)$Lon, sample_data(stem_otu)$Lat)
root.geo = cbind(sample_data(root_otu)$Lon, sample_data(root_otu)$Lat)
air.geo  = cbind(sample_data(air_otu)$Lon,  sample_data(air_otu)$Lat)
leaf.geo = cbind(sample_data(leaf_otu)$Lon, sample_data(leaf_otu)$Lat)
soil.geo = cbind(sample_data(soil_otu)$Lon, sample_data(soil_otu)$Lat)
litt.geo = cbind(sample_data(litt_otu)$Lon, sample_data(litt_otu)$Lat)
axil.geo = cbind(sample_data(axil_otu)$Lon, sample_data(axil_otu)$Lat) 
peti.geo = cbind(sample_data(peti_otu)$Lon, sample_data(peti_otu)$Lat)

# Generate a dissimilarity matrix for grouped (all plant parts) and by plant part 
#   bacterial geographic distance data 
all.geodist  = earth.dist(all.geo)
stem.geodist = earth.dist(stem.geo)
root.geodist = earth.dist(root.geo)
air.geodist  = earth.dist(air.geo)
leaf.geodist = earth.dist(leaf.geo)
soil.geodist = earth.dist(soil.geo)
litt.geodist = earth.dist(litt.geo)
axil.geodist = earth.dist(axil.geo)
peti.geodist = earth.dist(peti.geo)

#Mantel test by SampleType
all_bacteria_mantel  = mantel(log( all.geodist+1),log( all_otu.dist), permutations= 999)
stem_bacteria_mantel = mantel(log(stem.geodist+1),log(stem_otu.dist), permutations= 999)
root_bacteria_mantel = mantel(log(root.geodist+1),log(root_otu.dist), permutations= 999)
air_bacteria_mantel  = mantel(log( air.geodist+1),log( air_otu.dist), permutations= 999)
leaf_bacteria_mantel = mantel(log(leaf.geodist+1),log(leaf_otu.dist), permutations= 999)
soil_bacteria_mantel = mantel(log(soil.geodist+1),log(soil_otu.dist), permutations= 999)
litt_bacteria_mantel = mantel(log(litt.geodist+1),log(litt_otu.dist), permutations= 999)
axil_bacteria_mantel = mantel(log(axil.geodist+1),log(axil_otu.dist), permutations= 999)
peti_bacteria_mantel = mantel(log(peti.geodist+1),log(peti_otu.dist), permutations= 999)



#############################################
################ Fungal #####################
#############################################

#Load RDS file of fungal reads
FPS_normal = readRDS("data/fungal_PS_normal.gz") #"fungal_PS_normal"

#square root transform (Hellinger)
FPS_normal.hell = transform_sample_counts(FPS_normal, function(x) x^0.5)

#Subset Fungal Data by SampleType
fstem_otu = subset_samples(FPS_normal.hell, SampleType=="Stem")
froot_otu = subset_samples(FPS_normal.hell, SampleType=="Root")
fair_otu  = subset_samples(FPS_normal.hell, SampleType=="Air")
fleaf_otu = subset_samples(FPS_normal.hell, SampleType=="Leaf")
fsoil_otu = subset_samples(FPS_normal.hell, SampleType=="Soil")
flitt_otu = subset_samples(FPS_normal.hell, SampleType=="Litter")
faxil_otu = subset_samples(FPS_normal.hell, SampleType=="Axil")
fpeti_otu = subset_samples(FPS_normal.hell, SampleType=="Petiole")


#Generate Distance Matrix
fall_otu.dist  = phyloseq::distance(FPS_normal.hell, "bray")
fstem_otu.dist = phyloseq::distance(fstem_otu, "bray")
froot_otu.dist = phyloseq::distance(froot_otu, "bray")
fair_otu.dist  = phyloseq::distance( fair_otu, "bray")
fleaf_otu.dist = phyloseq::distance(fleaf_otu, "bray")
fsoil_otu.dist = phyloseq::distance(fsoil_otu, "bray")
flitt_otu.dist = phyloseq::distance(flitt_otu, "bray")
faxil_otu.dist = phyloseq::distance(faxil_otu, "bray")
fpeti_otu.dist = phyloseq::distance(fpeti_otu, "bray")


#Gather fungal spatial data from the phyloseq object
fall.geo  = cbind(sample_data(FPS_normal.hell)$Lon, sample_data(FPS_normal.hell)$Lat)
fstem.geo = cbind(sample_data(fstem_otu)$Lon, sample_data(fstem_otu)$Lat)
froot.geo = cbind(sample_data(froot_otu)$Lon, sample_data(froot_otu)$Lat)
fair.geo  = cbind(sample_data( fair_otu)$Lon, sample_data( fair_otu)$Lat)
fleaf.geo = cbind(sample_data(fleaf_otu)$Lon, sample_data(fleaf_otu)$Lat)
fsoil.geo = cbind(sample_data(fsoil_otu)$Lon, sample_data(fsoil_otu)$Lat)
flitt.geo = cbind(sample_data(flitt_otu)$Lon, sample_data(flitt_otu)$Lat)
faxil.geo = cbind(sample_data(faxil_otu)$Lon, sample_data(faxil_otu)$Lat)
fpeti.geo = cbind(sample_data(fpeti_otu)$Lon, sample_data(fpeti_otu)$Lat)


#Create correct sized fungal geographic distance dissimilarity matrix
fall.geodist  = earth.dist( fall.geo)
fstem.geodist = earth.dist(fstem.geo)
froot.geodist = earth.dist(froot.geo)
fair.geodist  = earth.dist( fair.geo)
fleaf.geodist = earth.dist(fleaf.geo)
fsoil.geodist = earth.dist(fsoil.geo)
flitt.geodist = earth.dist(flitt.geo)
faxil.geodist = earth.dist(faxil.geo)
fpeti.geodist = earth.dist(fpeti.geo)


#Mantel test by SampleType
fall_mantel  = mantel(log( fall.geodist+1),log( fall_otu.dist), permutations= 999)
fstem_mantel = mantel(log(fstem.geodist+1),log(fstem_otu.dist), permutations= 999)
froot_mantel = mantel(log(froot.geodist+1),log(froot_otu.dist), permutations= 999)
fair_mantel  = mantel(log( fair.geodist+1),log( fair_otu.dist), permutations= 999)
fleaf_mantel = mantel(log(fleaf.geodist+1),log(fleaf_otu.dist), permutations= 999)
fsoil_mantel = mantel(log(fsoil.geodist+1),log(fsoil_otu.dist), permutations= 999)
flitt_mantel = mantel(log(flitt.geodist+1),log(flitt_otu.dist), permutations= 999)
faxil_mantel = mantel(log(faxil.geodist+1),log(faxil_otu.dist), permutations= 999)
fpeti_mantel = mantel(log(fpeti.geodist+1),log(fpeti_otu.dist), permutations= 999)


###################################################
# Save Mantel Correlation Statistics and P Values #
###################################################

#Call plant part results
mant_cor = as.data.frame(c(stem_bacteria_mantel$statistic, root_bacteria_mantel$statistic, 
    air_bacteria_mantel$statistic, leaf_bacteria_mantel$statistic, 
    soil_bacteria_mantel$statistic, litt_bacteria_mantel$statistic, 
    axil_bacteria_mantel$statistic, peti_bacteria_mantel$statistic, 
    fstem_mantel$statistic, froot_mantel$statistic, fair_mantel$statistic, 
    fleaf_mantel$statistic, fsoil_mantel$statistic, flitt_mantel$statistic, 
    faxil_mantel$statistic, fpeti_mantel$statistic))

#Call plant part p-values
p_val = c(stem_bacteria_mantel$signif, root_bacteria_mantel$signif, 
    air_bacteria_mantel$signif, leaf_bacteria_mantel$signif, 
    soil_bacteria_mantel$signif, litt_bacteria_mantel$signif, 
    axil_bacteria_mantel$signif, peti_bacteria_mantel$signif, 
    fstem_mantel$signif, froot_mantel$signif, fair_mantel$signif, 
    fleaf_mantel$signif, fsoil_mantel$signif, flitt_mantel$signif, 
    faxil_mantel$signif, fpeti_mantel$signif)

#Create results data frame with mantel R and p-values for each plant part
results = data.frame( c( rep('Bacterial',8), rep('Fungal',8) ), 
    c( rep( c('Stem', 'Root', 'Air', 'Leaf', 'Soil', 'Litter', 'Axil', 'Petiole'),2 ) ),
    mant_cor, p_val
    )

#Name the columns intelligible names
colnames(results) = c('Type', 'Part', 'Mant_cor', 'p_val')

#Adjust p-values for repeated tests on the within plant part groups for bacteria
bact.cor.p = p.adjust(results[results$Type == 'Bacterial',4], method = 'bonferroni')
bact.cor.p = data.frame(cor.p = bact.cor.p, Type = rep('Bacterial',8), 
    Part = as.character(unique(results$Part)) )

#Adjust p-values for repeated tests on the within plant part groups for fungi
fung.cor.p = p.adjust(results[results$Type == 'Fungal',4], method = 'bonferroni')
fung.cor.p = data.frame(cor.p = fung.cor.p, Type = rep('Fungal',8), 
    Part = as.character(unique(results$Part)) )

#Join corrected p-values to for fungus and bacteria together
cor.p = rbind(bact.cor.p, fung.cor.p)


#Join corrected p-values to results data frame
results = merge(results, cor.p)


#Add overall test (all parts together) to data frame
results2 = rbind(results, data.frame(Type = c('Bacterial', 'Fungal'), Part = rep('All',2),
    Mant_cor = c(all_bacteria_mantel$statistic, fall_mantel$statistic), 
    p_val = c(all_bacteria_mantel$signif, fall_mantel$signif),
    cor.p = rep('NA', 2))
    )

#Save all result data as a csv
write.csv(results2, file = 'output/Mantel_Plant_Part_Results_Hell_Bray.csv')
knitr::kable(results2, digits = c(0, 0, 3, 3, 3, 3))
Type Part Mant_cor p_val cor.p
Bacterial Air 0.174 0.277 1
Bacterial Axil 0.075 0.334 1
Bacterial Leaf 0.213 0.122 0.976
Bacterial Litter 0.280 0.121 0.968
Bacterial Petiole 0.187 0.179 1
Bacterial Root 0.444 0.001 0.008
Bacterial Soil 0.044 0.384 1
Bacterial Stem 0.245 0.179 1
Fungal Air 0.613 0.001 0.008
Fungal Axil 0.225 0.186 1
Fungal Leaf 0.114 0.230 1
Fungal Litter 0.598 0.003 0.024
Fungal Petiole 0.670 0.001 0.008
Fungal Root 0.700 0.001 0.008
Fungal Soil 0.674 0.007 0.056
Fungal Stem 0.555 0.001 0.008
Bacterial All 0.116 0.004 NA
Fungal All 0.364 0.001 NA

PERMANOVA

PERMANOVA for BACTERIA

##PERMANOVA Bacteria
PS_normal = readRDS("data/PS_normal.gz")
PS_normal=transform_sample_counts(PS_normal, function(x) x^0.5)

PS_normal = subset_taxa(PS_normal, c(
                          Order!="Chloroplast", 
                          Kingdom!="Archaea",
                          Family!="Mitochondria",
                          Kingdom!="Unknown",
                          Kingdom!="Eukaryota"))


#Use base boxplot functin to look at sequencing depth
#boxplot(sample_sums(PS_normal)~sample_data(PS_normal)$FieldSite)
#boxplot(sample_sums(PS_normal)~sample_data(PS_normal)$SampleType)

#Bacterial ALL
Site <- as.factor(PS_normal@sam_data$FieldSite)
Sample_Type <- PS_normal@sam_data$SampleType
bact_ps <- phyloseq::distance(PS_normal, method="bray") #Normalized bact data - phyloseq-class
bact_meta_all = as(sample_data(PS_normal), "data.frame") #Extract sample data and coerce to df
bact_all = adonis(bact_ps ~ Site + Sample_Type, by = "margin", data = bact_meta_all, permutations = 10000)
bact_all
## 
## Call:
## adonis(formula = bact_ps ~ Site + Sample_Type, data = bact_meta_all,      permutations = 10000, by = "margin") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## Site         8    4.0127 0.50159  2.0569 0.15102 9.999e-05 ***
## Sample_Type  7    8.9018 1.27169  5.2150 0.33503 9.999e-05 ***
## Residuals   56   13.6557 0.24385         0.51395              
## Total       71   26.5702                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMANOVA for Fungi

##PERMANOVA Fungi
FPS_normal = readRDS("data/fungal_PS_normal.gz")

fun_ps = transform_sample_counts(FPS_normal, function(x) x^0.5)
site <- as.factor(FPS_normal@sam_data$FieldSite)
SampleType <- FPS_normal@sam_data$SampleType
fun_ps <- phyloseq::distance(fun_ps, method="bray") #Normalized bact data - phyloseq-class
fun_meta_all = as(sample_data(FPS_normal), "data.frame") # Extract sample data and coerce to df
fun_all = adonis(fun_ps ~ site + SampleType, by = "margin", data = fun_meta_all, permutations = 10000)
fun_all
## 
## Call:
## adonis(formula = fun_ps ~ site + SampleType, data = fun_meta_all,      permutations = 10000, by = "margin") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## site        8    6.5486 0.81858  3.2380 0.25541 9.999e-05 ***
## SampleType  7    4.9343 0.70491  2.7884 0.19245 9.999e-05 ***
## Residuals  56   14.1570 0.25280         0.55215              
## Total      71   25.6399                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NMDS

Non-metric multidimensional scaling (NMDS) using Bray-Curtis ordination for Hellinger transformed bacterial and fungal communities with elipses and points for each plant-part.

#Code by Maria and Helen
#For NMDS plots and heat maps by Sample Type

#read in normalized bacterial data
PS_normal<-readRDS("data/PS_normal.gz")
PS_normal<- subset_taxa(PS_normal, c(Kingdom != "unknown", Kingdom != "Eukaryota"))

#read in bacteria phyloseq data
physeq1 <- readRDS("data/physeq1.gz") 
physeq1<- subset_taxa(physeq1, c(Kingdom != "unknown", Kingdom != "Eukaryota"))

#read in rarified phyloseq data
rarPhySeq <- readRDS("data/rarPhySeq.gz") 
rarPhySeq<- subset_taxa(rarPhySeq, c(Kingdom != "unknown", Kingdom != "Eukaryota"))

#read in normalized fungal taxa
fungal_norm<-readRDS("data/fungal_PS_normal.gz")

#read in fungal phyloseq data
funphys <- readRDS("data/fungal_physeq1.gz") 

#read in rarified phyloseq data
rarfun <- readRDS("data/fungal_rarPhySeq.gz")

####################################
##########NMDS PLOTS##############
#hellinger transform bacteria
phyobj.hell=transform_sample_counts(PS_normal, function(x) x^0.5)

#hellinger transform fungus
fungal.hell<-transform_sample_counts(fungal_norm,  function(x) x^0.5)

##Organize sample type by levels of ground to air in all datasets
sample_data(PS_normal)$SampleType = factor(sample_data(PS_normal)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

sample_data(physeq1)$SampleType = factor(sample_data(physeq1)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

sample_data(rarPhySeq)$SampleType = factor(sample_data(rarPhySeq)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

##Organize sample type by levels of ground to air in all datasets for fungus
sample_data(fungal_norm)$SampleType = factor(sample_data(fungal_norm)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

sample_data(funphys)$SampleType = factor(sample_data(funphys)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

sample_data(rarfun)$SampleType = factor(sample_data(rarfun)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

#bacteria ordination
phyobj.ord<-ordinate(phyobj.hell, "NMDS", "bray")

#fungus ordination
funhell.ord<-ordinate(fungal.hell, "NMDS", "bray")
 
############## Black and White

#bacteria NMDS by sample type
bac.NMDS<-plot_ordination(phyobj.hell, phyobj.ord, type="samples", shape="SampleType", title="Bacteria") +
  theme(legend.position="none")+
  scale_shape_manual(values=c(0, 21, 2, 3, 4, 17, 15, 8),
                     name="Microhabitats")+
  #stat_ellipse(type="t", linetype=1, col="gray60", lwd=0.4) +
  coord_equal(ratio=2.68) +
  annotate(geom="text", x=-1.8, y=1.1, label=expression(bold("a")), size=3)+
  theme_bw()


#fungal NMDS by sample type
fung.NMDS<-plot_ordination(fungal.hell, funhell.ord, type="samples", shape="SampleType",  title="Fungi") +
  ylab("")+
  theme(legend.position="none")+
  scale_shape_manual(values=c(0, 21, 2, 3, 4, 17, 15, 8),
                     name="Microhabitats")+
  #stat_ellipse(type="t", linetype=1, col="gray60", lwd=0.4) +
  coord_equal(ratio=2) +
  annotate(geom="text", x=-1.2, y=1.1, label=expression(bold("b")), size=3)+
  theme_bw()

NMDS.legend <- get_legend(
  # create some space to the left of the legend
  fung.NMDS + theme(legend.box.margin = margin(0, 0, 0, 12)))

NMDS.plots<- plot_grid(bac.NMDS + theme(legend.position = "none"), 
                       fung.NMDS + theme(legend.position = "none"),
                       ncol=2,nrow=1)

plot_grid(NMDS.plots, NMDS.legend, rel_widths = c(3, 1)) # legend  1/3 size as first obj.

dev.copy(pdf, "figures/Fig2. NMDS.BW.pdf", width = 6, height = 4)

Heatmaps

####HEAT MAP DATA####

##Helen Sung, edited by Jared Bernard

###Load phyloseq files###
bacteria <- readRDS("data/PS_normal.gz") ## normalized phyloseq data
bacteria <- subset_taxa(bacteria, Kingdom != "unknown") #Remove unclassified kingdoms
bacteria <- subset_taxa(bacteria, Kingdom != "Eukaryota")   #Remove eukaryotes
bacteria <- subset_taxa(bacteria, Kingdom=="Bacteria")

fungi <- readRDS("data/fungal_PS_normal.gz") ## normalized phyloseq data
fungi <- subset_taxa(fungi, Kingdom=="Fungi")

###Fix Data to properly format for Heatmaps###
##Substitute all of the FieldSite 10 into 9 in all datasets 
gsub("10", "9", (sample_data(bacteria)$FieldSite))
gsub("10", "9", (sample_data(fungi)$FieldSite))

##Replace old fieldsite names with new fieldsite names in all datasets
sample_data(bacteria)$FieldSite 
sample_data(fungi)$FieldSite 

##Organize sample type by levels of ground to air in all datasets
sample_data(bacteria)$SampleType = factor(sample_data(bacteria)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))
sample_data(fungi)$SampleType = factor(sample_data(fungi)$SampleType, levels = c('Soil', 'Root', 'Stem', 'Axil', 'Petiole', 'Leaf','Litter','Air'))

###Subset Data###
# NOTE: Can put any physeq dataset (i.e. PS_normal, physeq1, rarPhySeq, phyobject.hell) into tax_glom(physeqobject, taxrank = "Rank you want to subset to") to subset
# tax_glom() gloms things together at the specified rank and it will get rid of anything that's not classified at that level

##Pool all ASVs by Class 
class_dat <- tax_glom(bacteria, taxrank="Class") ## subsetting for normalized data
funclass <- tax_glom(fungi, taxrank="Class")

###Plot Heatmap###
# NOTE: Can put any Subsetted physeq dataset (i.e. phy_dat, ord_dat, etc.) or non-subsetted physeq dataset (i.e. physeq1, PS_normal, etc.) into plot_heatmap(physeqobject, ...)
# NOTE: If using a subsetted physeq dataset, make sure set the taxa.label = taxrank for whichever taxonomic rank you make heatmap for

#Plot bacterial heatmap clustering by NMDS and Bray curtis distance, labeling sample type and supressing Taxon labels
b1 <- plot_heatmap(class_dat, "NMDS", "bray", sample.label = "SampleType", sample.order= "FieldSite", taxa.label="Class", low = "#3CBB75FF", high = "#440154FF", na.value= "white")
#Organize by field sites and grouped by sample type
d1 <- b1+ facet_grid(~SampleType, scales= "free_x", switch = "x")
#Plot Heatmap with tick marks and FieldSite labels removed
d2 <- d1 + theme(
  axis.text.x = element_blank(),
  axis.ticks = element_blank()
  ) 
plot(d2)

##Save Heatmap as tiff file
tiff("figures/FigS1. 16S_bact_heatmap.tiff", units="in", width=12, height=12, res=300) ## Change "ITS_order_heatmap.tiff" to whatever title for the file you choose
plot(d2)
dev.off()

####################

#Plot fungal heatmap clustering by NMDS and Bray curtis distance, labeling sample type and supressing Taxon labels
p1 <- plot_heatmap(funclass, "NMDS", "bray", sample.label = "SampleType", sample.order= "FieldSite", taxa.label="Class", low = "#FDE725FF", high = "#29AF7FFF", na.value= "white")
#Organize by field sites and grouped by sample type
q1 <- p1+ facet_grid(~SampleType, scales= "free_x", switch = "x")
##Plot Heatmap with tick marks and FieldSite labels removed
q2 <- q1 + theme(
  axis.text.x = element_blank(),
  axis.ticks = element_blank()
  ) 
plot(q2)

##Save Heatmap as tiff file
tiff("figures/FigS2. ITS_fungal_heatmap.tiff", units="in", width=12, height=12, res=300)
plot(q2)
dev.off()